/*
 * Created by edhong 
 * Mar 2008
 *
 * 
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License
 * as published by the Free Software Foundation; either version 2
 * of the License, or (at your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 *
 */


#include <math.h>
#include "libQccPack.h"
#define USG_STRING "[-bp %d:bitplanes] [-t %d:firstthreshold] %s:sbpfile_input %s:sbpfile_output"


int NumLevels ;

QccWAVSubbandPyramid sbpin, sbpout;

int numcols, numrows;
int temp;
double thresh;

int main(int argc, char *argv[])
{
  
  QccInit(argc, argv);

  thresh = 2048.0;
  NumLevels = 5;

  QccWAVSubbandPyramidInitialize(&sbpout);
  QccWAVSubbandPyramidInitialize(&sbpin);

  if (QccParseParameters(argc, argv,
                         USG_STRING,
                         &NumLevels,
			 &temp,
                         sbpin.filename,
                         sbpout.filename))
    QccErrorExit();

  thresh = temp;
      
  if (NumLevels < 1)
    {
      QccErrorAddMessage("%s: Number of bitplanes must be 1 or greater",
                         argv[0]);
      QccErrorExit();
    }

  if (QccWAVSubbandPyramidRead(&sbpin))
    {
      QccErrorAddMessage("%s: Error calling QccWAVSubbandPyramidRead()",
                         argv[0]);
      QccErrorExit();
    }

  printf("truncating %s with %d bitplanes, threshold of %f, output to %s\n",
	 sbpin.filename, NumLevels, thresh, sbpout.filename);

  numcols = sbpin.num_cols;
  numrows = sbpin.num_rows;

  int i,j;
  for (i=0; i < numrows; i++) 
    for  (j = 0; j < numcols; j++) {
      double currcoeff = sbpin.matrix[i][j];
      int sign;
      if (currcoeff <0) {
	sign = -1;
	currcoeff = -currcoeff;
      }
      else 
	sign = 1;

      if (currcoeff >= 2.0* thresh) {
	printf("Error: threshold too small for coeff (%d,%d) with magnitude %f\n", 
	       i,j, currcoeff);
      }
      double factor = thresh/pow(2,NumLevels-1);

      currcoeff = (floor((currcoeff/ factor)))*factor;      
      
      sbpin.matrix[i][j] = sign * currcoeff;
    }

  QccStringCopy(sbpin.filename, sbpout.filename);

  
  if (QccWAVSubbandPyramidWrite(&sbpin))
    {
      QccErrorAddMessage("%s: Error calling QccWAVSubbandPyramidWrite()",
                         argv[0]);
      QccErrorExit();
    }


  QccExit;
}
